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As a first step toward a fully two-dimensional asymptotic theory for the bifurcation of soli- 
tons from infinitesimal continuous waves, an analytical theory is presented for line solitons, whose 
envelope varies only along one direction, in general two-dimensional periodic potentials. For this 
two-dimensional problem, it is no longer viable to rely on a certain recurrence relation for going 
beyond all orders of the usual multi-scale perturbation expansion, a key step of the exponential 
asymptotics procedure previously used for solitons in one-dimensional problems. Instead, we pro- 
pose a more direct treatment which not only overcomes the recurrence-relation limitation, but also 
simplifies the exponential asymptotics process. Using this modified technique, we show that line 
solitons with any rational line slopes bifurcate out from every Bloch-band edge; and for each ra- 
tional slope, two line-soliton families exist. Furthermore, line solitons can bifurcate from interior 
points of Bloch bands as well, but such line solitons exist only for a couple of special line angles 
due to resonance with the Bloch bands. In addition, we show that a countable set of multi-line- 
soliton bound states can be constructed analytically. The analytical predictions are compared with 
numerical results for both symmetric and asymmetric potentials, and good agreement is obtained. 



1 Introduction 

Nonlinear wave propagation in periodic media is currently a subject of intensive research in optics 
and applied mathematics, with diverse applications ranging from nonlinear photonics [H[2] to Bose- 
Einstein condensates (3J H] . While many types of soliton structures in periodic media have been 
reported theoretically and experimentally (see [5] for a review), prior analytical work is limited to 
special (symmetric) periodic potentials P, [7] , or to one spatial dimension [U EJ [TO] . 

In this paper, we study line solitons in general two-dimensional periodic media, employing 
the nonlinear Schrodinger (NLS) equation with a spatially periodic potential as our model. Such 
solitons have been reported theoretically and experimentally in certain periodic potentials before 

IT!?! [TTH [T4"l I15| . and it is known that they can bifurcate from edges of Bloch bands or from 
certain high-symmetry points inside Bloch bands. However, systematic analytical construction of 
line solitons is still lacking, and a number of open questions remain. It is not clear, in particular, 
if the inclination angles and positions of these line solitons relative to the underlying potential can 
be arbitrary or not. 

*e-mail: jyang@math.uvm.edu 
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To address these issues, in this article we analytically examine small-amplitude line solitons 
which bifurcate out from infinitesimal Bloch modes. Such solitons are in the form of slowly varying 
Bloch-wave packets, with the envelope of the packet being uniform along the line direction. Using 
standard asymptotic methods, the packet envelope can be readily shown to satisfy the familiar one- 
dimensional NLS equation in the absence of a potential and thus have a sech-shape. However, the 
position of the envelope relative to the potential is harder to determine because it hinges on effects 
that are exponentially small in the soliton amplitude. For this purpose, techniques of exponential 
asymptotics must be used. 

In one-dimensional nonlinear wave systems, an exponential asymptotics method for analyzing 
low-amplitude solitons that comprise a periodic carrier modulated by an envelope, has been devel- 
oped in the past fifteen years [HI El \TU[ IT6] . This method focuses on the Fourier transform with 
respect to the slow spatial variable of the envelope function, motivated by the fact that the solitary- 
wave tails in the physical domain are controlled by the pole singularities near the real axis of the 
wavenumber space. The residues of these poles, which are exponentially small, are then calculated 
by matched asymptotics near the poles and away from the poles. The solution away from the 
poles, in particular, is expressed as a Taylor series of the wavenumber around zero, multiplied by 
a factor which is exponentially small at the poles, and the coefficients in this series are computed 
through a recurrence relation. Upon inverting the Fourier transform, these poles of exponentially 
small strength give rise to growing tails of exponentially small amplitudes in the physical solution. 
It turns out that these growing tails vanish only at two positions of the envelope relative to the 
periodic carrier, which in turn reveals that only two soliton families are admitted. 

In this paper, by investigating line solitons, we present the first step in developing a fully two- 
dimensional asymptotic theory for the bifurcation of solitons from infinitesimal continuous waves. 
For line solitons, the governing equation for the carrier envelope is still one-dimensional, hence the 
corresponding Fourier transform depends on only one complex variable in the wavenumber domain. 
While this simplifies matters, two new obstacles arise. The first is that, in a two-dimensional 
problem, the Fourier transform of the solution contains pole singularities not only near the real 
axis (which we call real poles), but also away from the real axis (which we call complex poles). 
Which of these poles are relevant for the determination of line solitons is unclear. On this issue, 
we will show that it is still the real poles which matter. 

The second and more serious obstacle is that, in two dimensions, the poles closest to the origin 
are often complex rather than real. As a consequence, utilizing a recurrence relation to determine 
the residues of real poles, a key step in one-dimensional problems, is no longer viable: the radius of 
convergence of the Taylor series (around the zero wavenumber) is limited by the nearest complex 
poles and thus cannot reach the real poles further away, which makes the matched asymptotics 
impossible to use. To overcome this obstacle, we shall give up the recurrence-relation approach and 
directly solve the Fourier-transformed equation along the real line of the wavenumber plane, up to 
the real poles. By doing so, complex poles become irrelevant and thus can be ignored. This modified 
approach not only turns out to be effective for two dimensions, but also simplifies the exponential 
asymptotics procedure: by working with the Fourier-transformed equation, the exponentially-small 
terms are factored out and one is left with a Volterra integral equation, which is in fact easier to 
solve numerically than the recurrence relation. 

We apply this modified exponential asymptotics technique to line solitons bifurcating from edges 
of Bloch bands in general two-dimensional periodic potentials. At low amplitudes, the solution is 
a Bloch-wave packet whose envelope varies only in the direction normal to the line soliton. By 
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using the Fourier-transform approach as outlined above, we analytically calculate growing tails of 
exponentially small amplitudes in the normal direction. The results show that for any rational 
slope, two line solitons relative to the potential exist, and the positions of their envelopes can be 
explicitly obtained. By matching the growing and decaying tails, we also construct an infinite 
family of multi-line-soliton bound states. These line solitons and multi-line-soliton bound states 
have been obtained numerically as well, in agreement with the theory. To quantitatively verify the 
analytical formula for the growing tails of exponentially-small amplitude, we use this tail formula to 
calculate the bifurcation of the zero eigenvalue associated with the linear stability of low-amplitude 
line solitons. This analytical eigenvalue formula is compared with the numerical eigenvalues, and 
good quantitative agreement is reached. Finally, we also show that line solitons can bifurcate from 
interior high-symmetry points of Bloch bands as well. But such line solitons exist only for very few 
(up to three) special line angles due to resonance with the Bloch bands. 



2 Line wavepackets at low amplitudes 

We consider the nonlinear Schrddinger equation in two spatial dimensions with a periodic (lattice) 
potential, 

m t + V 2 ^ -V(x,y)^ + a\^\ 2 ^ = 0, (1) 

where V 2 = d xx + d yy , the potential V(x,y) is periodic in both x and y, and a = ±1 is the sign 
of nonlinearity. This equation governs the nonlinear propagation of light in photonic lattices as 
well as the evolution of Bose-Einstein condensates in optical lattices [H [2J [3J [T7]. In the latter 
community, this equation is referred to as the Gross— Pitaevskii equation. 

In this article, we shall take the minimal periods of the potential V(x,y) to be the same in 
both x and y. This matching periodicity in x and y, while not necessary for our analysis, allows 
for simplification of the algebra. For a periodic potential with matching periodicity, rotation by an 
angle 0, with tan(#) rational, also yields a periodic potential, but the periods of the new potential 
would change in general. For instance, for the potential V(x, y) = sin 2 x + sin 2 y with matching 
periods it, a 45°-rotation would yield a potential with matching periods y2ir. To remove this 
rotational freedom and the changing periods of the same potential, we align the potential in such a 
way that its periods are minimal. We call the potential in such orientation as the minimal-period- 
orientation potential, and this potential is used throughout the paper. Without loss of generality, 
we also normalize the matching periods of this minimal-period-orientation potential as it. 

We search for stationary solutions to Eq. ([!]) of the form 

V(x,y,t) = TP(x,y)e- i »\ (2) 

where fi is the propagation constant and tjj is a real-valued amplitude function that satisfies the 
equation 

V 2 4> + (fi-V(x,y))4; + ai; 3 = 0. (3) 

In this article, we shall consider line-soliton solutions that are bounded (non-decaying) along a 
certain line in the (x, y) plane but decay to zero along the direction orthogonal to this line. These 
solutions may bifurcate out from edges of Bloch bands, or from certain high-symmetry points inside 
Bloch bands |12} fT3l [T5] . Our analysis will be developed first for line solitons bifurcating from 



edges of Bloch bands. Afterwards (in section 11) we shall consider line solitons bifurcating from 
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interior points of the Bloch bands, which do exist but under more restrictive conditions due to 
resonance with Bloch bands. 

Near edges of Bloch bands, solitons are low-amplitude wavepackets comprising the underlying 
Bloch- wave carrier modulated by a slowly varying envelope, and they can be studied using multiple- 
scale perturbation theory. For line solitons, the envelope only varies in one direction and the 
associated 'slow' spatial variable will be 

W = ew, w = x sin 9 — y cos 9, (4) 

where 9 is the angle the line soliton makes with the x axis, and y^/i — fio\ = e <C 1 is the deviation 
from the edge point /j,q. Plugging this into equation ^ yields 

L V> + edwL^ + e 2 d 2 v ip + atp 3 + {n~ Ho)^ = 0, (5) 

where 

[sin 6*, - cos#] . (6) 
The perturbation series solution to Eq. ^ is readily found to be 

il>(x, y, W) = eA(W)b(x, y) + e 2 A'(W)u(x, y) + 0(e 3 ), (7) 

where the functions of the fast variables satisfy 

L b(x,y) = 0, (8) 
L v(x,y) = -Lib(x,y). (9) 

Eq. ^ implies that b(x, y) is a Bloch mode at the edge point ^ = hq. For simplicity of the analysis, 
we shall make the following assumption. 

Assumption 1 The Bloch mode at the band edge /U = fio is unique (up to a multiplication 
constant). 

This assumption is for the purpose of avoiding nonlinear interactions between different Bloch 
modes at the same band-edge point, and it holds for many band edges. We also normalize b(x, y) 
so that max|6| = 1. 

The inhomogeneous equation ^ is solvable, since the right-hand side —L\b is clearly orthogonal 
to the homogeneous solution b{x, y) so the Fredholm condition is satisfied. To avoid ambiguity of 
the homogeneous term in u, we require {v, b) = 0, where the inner product is defined as 

(/, g) = r r f* ( x ' y^ x > y) dxd ^' ( 10 ) 

Jo Jo 

with the superscript '*' representing complex conjugation. 

From the solvability condition at order e 3 , the envelope function A(W) in expansion ([7]) satisfies 

DA" + r]A + aaA 3 = 0, (11) 

with 

(Lxv + b^) (b 3 ,b) 

D= (b,b) ' v = sm(v-vo), « = W - (12) 
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When sgn(cr) = sgn(D) = — sgn(r/), the solution of the envelope equation (11) is 



W - Wn 

A(W) = asech ° , (13) 

y 

where W = Wq is the center position of the line envelope, and 

a = <s/2/a, = y/\D\. (14) 

The asymptotic expansion Q may be carried to all powers of e, with all terms being localized in 
W (i.e., localized along the normal direction of the line solution), suggesting that line solitons exist 
for all choices of Wq. However, based on previous experience and numerical computations, this is 
not the case. Specifically, it turns out that for rational tan#, only two line-soliton families exist 
(not counting their periodic replications in the lattice). These two soliton families are illustrated 
in Fig. [I] for a = 1, e = 0.25 and tan# = 1, 2 in the specific lattice 

V(x, y) = 6 (sin 2 x + sin 2 y) , (15) 

where they bifurcate out from the lowest Bloch-band edge /io = 4.1264. Solitons in the upper 
panels have Wq = and are called onsite solitons, while those in the lower panels have Wq = eir/2 
and are called offsite solitons. 

This discrepancy between the prediction of the multi-scale expansion Q and true solutions was 
first studied in [18] for a one-dimensional problem, and it was suggested that the issue could be 
resolved by requiring that a Melnikov-integral condition be satisfied by the leading-order approxi- 
mation of ((7]). While this constraint happens to specify Wq correctly for symmetric potentials, in 
a general lattice all terms in the perturbation expansion ([7]) make contributions to the Melnikov 
integral at the same order of e, and the calculations involved quickly get out of hand. To handle 
this difficulty, an exponential-asymptotics approach was developed for one-dimensional problems in 
[SI [10] , following an earlier treatment of a similar problem in the fifth-order KdV equation |16j . In 
the next sections, we shall adapt this exponential-asymptotics method to the study of line solitons. 
For this two-dimensional problem, certain key steps in the previous exponential-asymptotics pro- 
cedure are no longer viable, and suitable modifications will be needed in order to overcome those 
obstacles. 



3 Solution in the wavenumber domain 

The failure of the multi-scale expansion ([7]) stems from the fact that, for generic values of the 
envelope position Wq, the exact solution ip(x,y;W) contains growing tails when W <C — 1 or 
W S> 1; but the amplitudes of these tails are exponentially small in e, and hence invisible in the 
expansion ([7]). In order to obtain true line solitons, we shall first calculate these exponentially small 
growing tails and then determine the envelope position Wq by insisting that these tails vanish. 

As before [HI \TU[ I16j . we shall work in the wavenumber domain. First, we take the Fourier 
transform of ip(x, y, W) with respect to the slow variable W, 

~ 1 f°° 

^{x,y,K) = — i;(x,y,W)e- lKW dW. (16) 

Under this transformation, the growing tails of exponentially small amplitude in ip(x, y, W) are 
transformed into poles of ip with exponentially small residues. 
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Taking the Fourier transform of the series solution ([7]) term-by-term yields 



a/3 



AWqK 



ip(x,V,K) = e—e x } 



scch 



7T/3K 



{b(x,y) + \eKv(x,y) + ...}. 



(17) 



This series is disordered at eK ~ 1. Thus, we introduce the 'slow' wavenumber k = eK and 
rearrange this series as 



where 



U(x,y,K;e) 



( e -iWW, sech ( ^] U{x,y,K;e) 



a(3 



2e 



{b(x,y) + \Kv{x,y) + . . .} , k < 1. 



(18) 



(19) 



We now derive the governing equation for U by taking the Fourier transform of equation ^ to 
arrive at 

Lq^+ikLi^ - K 2 tjj + aip 3 + e 2 r] ; tjj = 0. (20) 



Then, by substituting the expression (18), we find that U satisfies 



L U + ikLxU + (rje 2 - n 2 )U + a cosh 



tt(3k 
~2e 



oo roo 



U{k-t)U(t - s)U(s) 



- x ,'- x ( , )S h ( ZEfcl) cosh M=4) cosh 



drds = 0. 



(21) 



4 Poles in the wavenumber plane 

We are concerned with pole singularities in U (x, y, k; e) which account for the growing tails in the 
physical space. Singularities of U are expected to occur near values of k = kq where the linear part 



of equation (21) is zero, i.e., 



L (p + iK Li(j) - nl4> = 0. (22) 

With a change of variables <f) = e~ 1K ° w (f>, where w is defined in Eq. equation (22) reduces to 

L O = O, (23) 

which has a single bounded solution <p = b(x,y), the Bloch mode at band edge hq. Thus, if we 



restrict ourselves for the moment to real values of Kq, bounded solutions to Eq. (22) are 

- iK0W b(x,y). 



(24) 



Since the spatial period of the solution <f>(x,y) should match that of the solution (19), <j)(x,y) and 



b(x, y) should have the same periodicity in (x, y). Then following the argument in jS], both kq cos(#) 
and Kosin(#) should be even integers. For this to occur, tan# must be a rational number, say 

tan(0)=p/?, (25) 
where p and q are relatively prime. Now to satisfy the periodicity condition, we get 

k = 2nvV + q 2 (26) 
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for any integer n. Thus poles of U(x, y, k; e) near the real axis of k are located near these kq values. 



To get the approximate locations of all poles in U(x, y, k; e), we repose Eq. (22) as an eigenvalue 
problem 



-1 



(27) 



where 



] T , with the superscript 'X" denoting vector transpose. It can be easily 
shown that all eigenvalues in (27) come in quadruples (kq, — kq, Kq, — Kq). This spectrum also has 



(l/i«o) 



periodicity 2y / p 2 + q 2 in the real direction of kq 



Numerically we solve the eigenvalue problem (27) using the Fourier collocation method to get 



rough estimates of the spectrum, followed by the Newton-conjugate-gradient method to calculate 
particular eigenvalues to high accuracy [5]. Examples of this spectrum are shown in Fig. |]for 
the lattice (15) at fiQ = 4.1264 (the edge of the semi-infinite gap) and tan# = 1,2. Notice that 



this spectrum contains not only the real eigenvalues given by Eq. (26), but also a large number of 



complex eigenvalues. In addition, when tan# ^ 0, the eigenvalues closest to the origin are complex 
eigenvalues rather than real eigenvalues. 
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Figure 2: Locations of singularities in the wavenumber domain for the lattice (15) with tan# = 1, 2 



The spectra in Fig. [2] raise serious questions on the applicability of the exponential-asymptotics 
method, as used earlier for solitons in one-dimensional lattices, to this two-dimensional problem. 
We recall that, in the one-dimensional case [3[T0l[T6], this spectrum contains only real eigenvalues; 
and it is those real eigenvalues which specify the envelope positions. In the present case, where 
both real and complex eigenvalues exist, which of those eigenvalues are responsible for selecting the 
envelope positions? 

To tackle this issue, we consider line solitons in a simpler stripe lattice, where the lattice varies in 
the x-direction only. The particulars of the analysis are covered in the appendix. For such a lattice, 
we find that the spectrum contains only complex eigenvalues and no nonzero real eigenvalues. Here, 
however, the envelope position of a line soliton can be arbitrary, as solutions with different envelope 
positions are equivalent to each other under a vertical translation. This indicates that the envelope 
position is not affected by complex eigenvalues, and thus we shall focus on poles of U(x, y, k; e) 
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near the real eigenvalues (26). In particular, we need to determine the residues of poles near the 
smallest real eigenvalues 

Ko = ±2^/p 2 + q 2 , (28) 

since those poles provide the dominant contributions to the solitary-wave tails. 

Even though only real eigenvalues affect the envelope position, the fact that the eigenvalues 
closest to the origin are complex rather than real, as shown in Fig. [2j poses a major obstacle to 
the calculation of the residues of poles near the real eigenvalues kq = ±2y / p 2 + q 2 by the previous 
exponential-asymptotics technique [HJ [TUJ [16] . This difficulty and its resolution will be elaborated 
in the next section. 



5 Solution away from the poles 



The basic idea for calculating the residues of poles in U(x,y,K;e) is to match the 'inner' solution 
near the poles to the 'outer' solution away from the poles. Since the poles of interest are near 
kq = ±2y / p 2 + q 2 , in this section we will determine the solution U(x, y, k; e) for real values of k in 
the interval — kq < k < kq but not close to ±Ko- 

When e — > 0, the main contribution of the double integral in (21) comes from the triangular 



region < s < k, < s < r when k>0otk<s<0, r<s<0 when k < 0. Over this region, 



cosh 



it (3 k 



cosh 



7T/3(k 



2e 



cosh 



7r/3(r 



2e 



cosh 



it (3 s 
~2e 



(29) 



while outside this region the same expression is exponentially small. Thus, to 0(e 2 ) the integral 



equation (21) reduces to the "outer" equation 



LaUW + iKLiUW-^uW 

+ 4a dr U w (x,y,K-r) ds C/ (0) (x, y, r - s)U {0 \x, y, s) = 0, 
J Jo 



(30) 



where U(x,y, k; e) = U^°\x,y, k) + 0(e 2 ). Here the main contribution to the error of the integral 



comes from the three corners of the triangular integration region in the integral of (30); thus 
this error is 0(e 2 ), a result that has also been checked numerically for randomly chosen analytic 
functions of U(k). 

In previous applications of exponential asymptotics, this outer integral equation was solved by 
expanding U^°\x,y, k) into a power series of k, which turns the outer integral equation into a 
recurrence relation for the coefficients of the power series [SI fTU| [TB] . This treatment was possible 
since the poles of interest were closest to the origin, hence the power series was convergent up to 
those poles. In the present two-dimensional problem, however, this treatment fails, because the 
power series for U^(x,y,K) has radius of convergence bounded by the distance to the nearest 
complex poles. Thus this power series cannot converge near the real poles of interest. As a result, 
the previous approach of relying on the recurrence relation needs to be abandoned. 



Instead, we propose to solve the outer integral equation (30) numerically. Since this is a Volterra 
integral equation, it can be easily tackled by explicit numerical methods. First we discretize k, 
K n = uAk, and write 

U ( -°\x,y, Kn ) = U n (x,y). (31) 
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Then we approximate the integral in (30) using the trapezoid rule. After the terms in the resulting 

(32) 



equation are rearranged, U n is found to satisfy a linear inhomogeneous equation 

[Lo + \K n L x -K 2 n + 2<jU 2 {Ak) 2 ] U n = -4aF n , 
where the inhomogeneous term F n is given by 

n-l 



m— 1 



I 

-^UqIu + U m I n -, 

m=l 

I m = U l U m-l, for 1 < m < n, 
1=0 
n-l 

In = J2 U i U n-l- 



(33a) 
(33b) 
(33c) 



i=i 



Since the inner sums I m for m < n do not change on further iterations, these only need to be 



computed once. Notice that the homogeneous operator in Eq. (32) is self-adjoint, thus this linear 



inhomogeneous equation can be solved by the preconditioned conjugate gradient method. The 



the equation ( 19 ) as 



initial conditions Uq and U\ cannot be derived from Eq. (32) itself, but they can be obtained from 

a/3 



U (x,y) 



a/3 



b(x,y), Ui{x,y) 



[b(x,y) + iAk v(x,y)] 



(34) 



The above numerical scheme for solving the outer integral equation (30) is explicit. Our nu- 



merical testing shows that its numerical error is O(Ak), thus it is first order accurate in Are. If 
one wishes for a higher-order numerical scheme, then instead of the trapezoidal rule, one can use a 
higher-order quadrature method (such as Simpson's rule) to approximate the integral in ( |30[ ). 

From the local analysis near the poles in the next section, it will transpire that U^°\x, y, k) has 
a fourth-order pole at Kq. Specifically, 



x,y,K) 



12C b(x,y) 



-IKQW 



as k 



5 {k-k Y 

where C is a complex constant. With a change of variables 

U(x,y,K) = (k- KofU^ix^,^, 



then 



U(x,y,K) 



12C 



b(x,y)e- iK ° l 



as k 



Kq. 



(35) 



(36) 



(37) 



Numerically we have confirmed the above outer-solution behavior near the poles. As an example, 
the numerical results for the lattice (15) with a = 1, tan# = 1 and /j,q = 4.1264 (edge of the semi- 
infinite gap) are shown in figure [3J For this line slope, kq = 2\/2. On the top left the contour plots 
of U(x,y, kq) are displayed. The corresponding analytic formula (37) with real C is also shown at 
the bottom for comparison. One can see that the two solutions match very well. On the right side 
of Fig. [3| the solution U (0, 0, k) is shown. When k — » kq, U(0, 0, k) — > 22X15, thus the fourth-order 
pole at k = acq is numerically confirmed. If tan# = 2, then we find that U(0, 0,k) — > 4.04 x 10 as 
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k — y kq = 2yE. Recalling the analytical formula (37) and Bloch-mode normalization (which boils 
down to 6(0,0) = 1 here), the constant C is then inferred as 



C = 9.19, for tan0 = 1, (38) 

C = 1.68xl0 4 , for tan 9 = 2. (39) 

In section [9] we are able to verify these values of C quantitatively, by comparing our predictions to 
numerical calculations of a certain linear-stability eigenvalue. 



Numerical U: Magnitude Phase 




-10 1 -10 1 K 

X X 



Figure 3: (Color online) Numerical solutions of the outer integral equation (30) for the lattice (15) 
with a = 1 and tan# = 1. Left panels: solutions at the singularity k = kq; upper row: the numerical 
solution, lower row: the analytical solution. Right panel: the solution versus k at x = y = 0. 



The approach taken here of directly solving the outer integral equation ( 30 ) by numerical meth- 
ods not only overcomes the inadequacy of the recurrence relation (for two-dimensional problems), 
but also has a number of additional advantages. Specifically, compared with the computation of 
recurrence relations in earlier works [TO] [T6] , the direct numerical solution of the outer integral 
equation, as was explained above, is actually easier. In addition, this direct computation gives di- 



rectly the pole strength (35) in the outer solution, while in the previous approach the pole strength 
was inferred indirectly from the recurrence solution. In view of these advantages, it is concluded 
that this direct solving of the outer integral equation also simplifies the exponential asymptotics 
procedure. 



6 Solution near the poles 

Based on experience from prior work [HI [10], [16] , the actual poles in the solution U(x,y,n;e) are 
expected to be O(e) away from the real axis. To determine the residues of these true poles, we 
need to analyze the solution U(x, y, k; e) near these poles. In this 'inner' region, the reduced outer 



equation (30) does not hold, and one needs to work with the full equation (21) instead. For line 
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solitons, it turns out that the analysis of the inner solution nearly replicates that in earlier works 
[HI HOj IS] , thus we shall keep this discussion brief. 

Focusing on the behavior of the solution U(x, y, k; e) in the inner region near the closest positive 
singularity, 

K = 2 V / p 2 + q 2 , (40) 
we introduce an inner variable, £ = (k — Ko)/e, that is k = kq + e£, with £ = 0(1). In this region 



we expand the solution to integral equation (21) as 

p — in w 



U(x,y,n;e) 



{$o(e)6(x, y) + ie&otiHx, v) + 0{e 2 )} . 



(41) 



Here the order of the solution near the pole U = 0(e 4 ) is chosen to match the large-£ behavior of 
the solution $ Q (£) = 0(£~ 4 ) (see Eq. $ty below). 



The dominant contribution from the double integral in equation (21) comes from the three 



regions: (i) r ~ 0, s ~ 0, (ii) r ~ k, s ~ 0, and (iii) r ~ k, s ~ k. In the first region, 



U(x, y,r — s) 
U(x,y,K- r) 



a/3 

U(x,y,s) « — b(x,y), 



-$o ( ; ) b(x,y), 



7t/3k / tt(3(k — r) 
cosh — — / cosh ■ 



(42) 
(43) 
(44) 

(45) 



the contribution from the double integral in the first region can be readily obtained. Contributions 
from the other two regions can be calculated in a similar manner, and they turn out to be the same 
as that from the first region. 

With the change of variables U(x, y, £) = e~ lKoW U(x, y, £), we then find that near the singularity 



2e / 2e 
Changing variables s = s/e, r = r/e and using the fact that 

sech(x — ?/)sech(y) dy = 2xcsch(x), 



the full integral equation (21) reduces to the inner equation 

3 

2e 2 



L U + etUXJ + e 2 {7] - f)U + ^aa 2 /3 2 6(x, yf 

/ tt/3o; 



X 



DC 



^/ 2 csch 



<I>o(C - w)dw = 0. 



(46) 



After substituting in the expansion (41), at 0(e 4 ) and 0(e 3 ) the above equation is automatically 
satisfied; and at 0(e~ 2 ) the equation governing 3>o(0 is found from the solvability condition as 



(1 + ^f)^ - 3 /3 2 y°° we ^/ 2 csch f 7 ^) $ (£ - ^du = 0. 



(47) 



This is a linear homogeneous Fredholm integral equation which has been solved before [HI [16] . Its 
analytic solution in the region |Imag(£)| > 1//3 is 

6/3 4 '• ±i0 ° 



MO 



i + p 2 e 



sin 2 s 



is e 



(48) 
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where the plus sign in the contour is for Imag(£) < — 1//3, the minus sign for Imag(£) > 



C 



2 cos 2 s 3s cos s 

: + -: 



sin 2 s 



(49) 



sins sins 

and C is a complex constant. Clearly this solution has simple-pole singularities at £ = ±i//3. Since 
the integral of (48) at these points is equal to — C/6, we see that 



(50) 



thus the pole has strength ±i/3 3 C/2 at £ = ±i//3. After substituting this back into equation (41) 
and changing variables to K = kq/e + ^, we find that U(x, y, K) has simple poles at K = Ko/e±i/ f3, 
and 



U(x,y,K) 



-1KQW 



b(x,y)- 



1 



K 



KO _l_ _1_ 



for if 



^0 , 1 



(51) 



Finally, from Eq. (18), we obtain the local behavior of tf)(x,y,K) near the simple poles K 

Ko/e±i/f3 as 



ij{x,y,K) ~ ^ e -T0Ko/2« e ±W>//9_?_^ r 6(x,y), K 

where 

wq = W /e. 

Moreover, from the symmetry of the Fourier transform 

V, K) = ^*{x, y, -K*) 



K0 l 



(52) 



(53) 



for real functions rp(x, y, W), we also deduce the local behavior of ip(x, y, K) near the simple poles 
K = -Ko/e±i/p. 

From the above local analysis, the residues of these poles are only determined up to a constant 
multiple (C is unknown at the moment). To determine C, we match the large-£ asymptotics of 
the above inner solution for U with the outer solution for U away from the singularities, in the 
matching region 1 <C |£| *C e _1 - To this end, we note that for |£| ^> 1 the main contribution of the 
integral in solution (48) comes from the region s ~ where 4>(s) ~ fC^ 3 - This yields 



1 1(1 1 

*o(fl~-5-^, for|£|»l. 



(54) 



Putting this back into equation (41), we find the following large-£ behavior for the inner solution, 
As discussed in section [5j this behavior matches to the outer solution with | k — kq | <C 1 . 
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7 Inversion of Fourier transform and true line solitons 



We now take the inverse Fourier transform of (16) 



if)(x,y,W) = J ip(x,y,K)e iKW dK (56) 
in order to determine the tail behaviors in the physical solution tp(x,y,W). Here ip(x,y,K) is 



given by Eq. (18). As explained in [8j, if we require this physical solution to decay upstream 
(w — > — oo), then the contour C in this inverse Fourier transform should be taken along the line 
Imag(K) = —1/(3 and pass below the poles K = ±Ko/e — i/ j3. It should also pass above the pole 



K = —i/(3 of the sech (irf3K/2) term in Eq. (18). Then when w ^ 1 (downstream), by completing 
the contour C with a large semicircle in the upper half plane, we pick up dominant contributions 
from the pole singularities at K = ±Ko/e — i/(3 and K = i/ j3. Collecting these pole contributions, 
the wave profile of the solution far downstream is then found to be 

■4> ~2eae- {w - Wo)/ Pb(x,y) 
4ir(3 3 C 



e 3 



-*0*o/ 2e s in{ Ko w - Q )e ( - W - Wo)/I3 b(x, y), w > 1/e, (57) 



where Co > and 0o are the amplitude and phase of the constant C (which is complex in general) , 
i.e., C = C e ie °. 



For this solution to be a line soliton, the growing term in (57) must vanish so sin (kowq — ©o 



0. Thus, there are two allowable locations for line solitons (relative to the lattice), 

wo = 6 /ko, (7r + e )/K . (58) 
Line solitons at these two locations are called onsite and offsite solitons respectively. For the 



particular lattice (15) and a = 1, these two line solitons near the edge of the semi-infinite gap with 
tan# = 1, 2 are displayed in Fig. [T] For these solitons, Oo = since C is real positive (see equations 
(38} and ([39])). 



The above results show that for any rational slope tan 9, two line solitons with envelope locations 
( 58 ) exist in a general two-dimensional lattice. What if the slope is irrational? Treating an irrational 



number as the limit of a rational number p/q with p, q — > oo, then kq = \/p 2 + q 2 — > oo, hence the 



growing tail downstream in (57) vanishes. This suggests that for an irrational slope, line solitons 
exist for arbitrary envelope positions wo. But we cannot numerically verify this conjecture since 
irrational numbers cannot be represented accurately on computers. 

Finally, if one wishes to obtain line wave packets ip(x,y, W) which decay for w — > +oo but 
contain a growing tail for w — 1, then the contour C in the inverse Fourier transform should be 
taken along the line Imag(K) = 1/(3 and pass above the poles K = ±Ko/e + i/f3 and below the pole 
K = i/(3. Then when w —1, by completing the contour C with a large semicircle in the lower 
half plane and picking up dominant pole contributions, the wave profile of the solution is found to 
be 

i(j ~2eae (w - WoW b(x,y) 

- 47r ^ C ° e -/W2e sin{KQWo _ e-oK^-^&Or, y), w < -1/e. (59) 

This 'flipped' wave solution will be useful when we construct multi-line-soliton bound states in the 
next section. 
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8 Construction of multi-line-soliton bound states 



The asymptotic tail formula (57) can be used not only to determine the locations of line solitons, but 
also to construct multi-line-soliton bound states. These bound states are analytically constructed 
by matching the downstream growing tail of a line wavepacket with the upstream decaying tail 
of another line wavepacket. This technique has been used for the construction of one-dimensional 
multi-soliton bound states before EH US] . Here we apply the same principle to the construction 
of multi-line soliton solutions in a general two-dimensional lattice potential. 

Consider the asymptotic expansion of two line wavepackets centered at w\ = W\ / e (left) and 
w 2 = Wilt (right) receptively. The left wavepacket decays for w — w\ <C — 1 and has a growing 
exponential tail for w — w\ ^> 1, and the right wavepacket decays for w — w 2 3> 1 and has a growing 
exponential tail for w — w 2 <C — 1. For line-soliton bound states, the decaying and growing tails 
of the two wavepackets must match in the region i»i < w < i«2. In this matching region, the 
left wavepacket's asymptotics is given by (57) with wq replaced by w\, and the right wavepacket's 
asymptotics is given by ( 59 ) with wq replaced by w 2 • Matching of these asymptotics results in the 
following system of equations 

2eae-( w - w ^b(x, y) = T fZ^V^°/ 2e sin (k w 2 - 6 ) e^ w - w ^b(x, y), 

2ea^ w - w ^b(x,y) = ±^^V^«o/2e sin ( Kom _ q q) ^W-W,)^ } 



where the =F comes from the possible it phase shift between the two wavepackets. After simplifica- 
tion, these matching conditions read 

sin (k wx - Go) = - sin (k w 2 - 6 ) = i^^e^/V^™^. (60) 

With a change of variables w\ = — {w\ — Qq/kq), w 2 = w 2 — Go/^o, the above matching conditions 
then become almost identical to the ones derived in [9] before. As has been explained there, this 
system of equations admits an infinite number of solutions for each fixed e > 0. Varying e, then 
infinite families of two-line soliton bound states are obtained. 

Note that for equations (60) to have a solution, the right-hand side must have magnitude less 
than one. Since this is not the case as e — > for any finite distance w 2 — w\ between the two 
line wavepackets, bifurcations of these bound states occur at finite amplitude away from the band 
edge. In figure |4j we show a particular family of bound states in the lattice (15) with a = 1 and 
tan# = 1. This family bifurcates at \i ~ 4.073 (near the edge [iq = 4.1264 of the semi-infinite gap). 
As predicted by the analysis of equation (60) in [9j, this solution family contains three connected 
branches, and their power curve is shown in Fig. [4] (upper left panel). Here the power is calculated 
over one period along the line-soliton direction w' = xcosQ + ysinO, which is ir\/p 2 + q 2 in the 
present case. Profiles of bound states at fj, = 4.0639 of the three power branches are displayed in 
Fig. [4](A,B,C). The bound state in the A panel comprises roughly two onsite line solitons, the one 
in the B panel comprises roughly an onsite and an offsite line solitons, and the bound state in the 
C panel comprises roughly two offsite line solitons. 

We would like to point out that the above construction of infinite families of line-soliton bound 
states was performed for general two-dimensional lattices. This contrasts the earlier work in [91 [10] 
where such construction was made only for symmetric one-dimensional lattices (where Go = 0). 
From the above calculations, it is now clear that the previous derivation of multi-soliton bound 
states in [9j [10] can be extended to general one-dimensional lattices, too. 
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Family of Bound States 



A 




9 Connection to zero-eigenvalue bifurcations 



In this section, we consider the linear stability eigenvalues of single-line solitons as obtained in 
section [7j Our interest lies in the pair of exponentially small eigenvalues of these solitons, which 
bifurcate out from the origin at the band edge; the associated eigenfunctions have zero wavenumber 
along the line-soliton direction (i.e., are periodic along the line direction with the period matching 
that of the line soliton). Such eigenvalues can be analytically calculated from the tail asymptotics 



(57) of line wavepackets which we have derived. Thus, by comparing the analytical formula of these 



eigenvalues with the numerically computed eigenvalues, we can quantitatively verify the formula in 



Eq. (57) for the exponentially small growing tails. We do caution, however, that we are ignoring 
other eigenvalues whose eigenfunctions have nonzero wavenumber along the line-soliton direction, 
and those eigenvalues are often more unstable [144 H5j . Thus the eigenvalue calculation in this 
section does not constitute a full stability analysis. 

The present calculation of zero-eigenvalue bifurcation closely parallels that in one-dimensional 
problems [191 E ttH]- Let ip s (x,y) = if)(x,y;wo s ) be a single-line soliton solution of equation Q 
with center at wo = wo s , which decays to zero as w — > =too and is periodic along the line direction, 
w' = x cos 9 + y sin 9, with period matching that of the Bloch wave. Perturbing this line soliton by 
normal modes 

* = e"*"* [V s + (v + p)e xt + (v* - ^V* 4 ] , (61) 
with v, tp <C 1, we obtain the linear-stability eigenvalue problem 

Codv = -X 2 v, (62) 

where 



C = V 2 + n - V(x, y) + a^ 2 s , A = V 2 + /x - V(x, y) + 3a^, 



and A is the stability eigenvalue. Since the bifurcated eigenvalue A is small near band edges, we 
expand the eigenfunction v into a perturbation series 

v = v + \ 2 V! + X 4 v 2 + ■■■■ (63) 



Inserting this expansion into Eq. (62), at 0(1) we get 

L L lVo = 0. (64) 

A solution to this equation is 

v = (dip/dw ) Wo=Wos ■ (65) 
Recalling the perturbation series expansion of if)(x, y; wq) in Eq. ([7]) as well as the envelope solution 



in Eq. (13), we find that 



e 2 a e(w-w 0s ), e(w - w Qs ) . . 
v r ^ sech p tanh p o{x,y). 

From this equation it is seen that the eigenfunction v is periodic along the line-soliton direction 
w' = x cos 9 + y sin 9, with period matching that of the Bloch wave (and the line soliton), so the 
net wavenumber of this eigenfunction along the line-soliton direction is zero. 
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From the large-u; asymptotics (57) of the solution ip(x,y;wo), we see that when w S> 1/e, fo 
contains a growing tail, 



v 



±e~ 3 4TTf3 3 K C Q e- 7T ^ /2e e e(w ~ w ^^b(x, y) . 



(66) 



Here the plus and minus sign correspond to the onsite and offsite line solitons, respectively. This 



growing tail must be balanced by the higher-order terms in the expansion ( 63 ) , thereby yielding an 



analytical formula for the eigenvalue A. Carrying out this program as in [51 110| [T9]. we find that 



A 2 = TC, 



ac 



(67) 



Thus this eigenvalue is stable for onsite line solitons and unstable for offsite ones, and its magnitude 
is exponentially small. 



A comparison of the analytical prediction (67) for the eigenvalue to numerical results for offsite 
line solitons is presented in Fig. [5] for the lattice (15) with a = 1 and tan# = 1, 2. It is seen that 
the analytical and numerical eigenvalues agree with each other very well. Notice that the analytical 
eigenvalue formula contains the constant Co, and the ratio of A 2 /[327TKo/3 4 e _7r ^ KO//2e /ae] approaches 
Co when e — > 0. We have plotted this ratio for numerically obtained eigenvalues in the right panels 
of Fig. [5] It is seen that as e — > 0, this ratio indeed approaches the C value obtained from equations 
(38]) and ([39]) in section [H] Thus our formula ([57]) for the exponentially-small growing tails in line 
wavepackets is fully verified quantitatively. 



10 Further examples with asymmetric potentials 

To help demonstrate the generality of our theory, we now consider the asymmetric lattice 

3 

V(x, y) = -(sin2x + sin4x + sin2y + sin4y), (68) 

which is shown in Fig. [6] (top left panel). For this lattice, the edge of the semi-infinite gap is 
/io = —0.6891. With soliton inclination tan# = —1 and nonlinearity coefficient a = 1, we find the 
two soliton families which bifurcate off this band edge (see Fig. [6] bottom panel). In this case 
numerical solution of the "outer" integral equation ( 30 ) yields a complex value for C = Coe 100 , 



C = 8.25, 9 = 2.66, (69) 
see Fig. k3] (top right panel). This 0q value in turn gives the center of the solitons, from formulae 



(58), as wo = 0.94 and 2.05 for onsite and offsite, respectively (up to periodic repetitions in the 
lattice) . Comparison of this prediction with the envelope locations of numerically obtained solitons 
shows good agreement (see Fig. [6] bottom panel). 



11 Line solitons bifurcated from interior points of Bloch bands 

In previous sections, we studied line solitons that bifurcate from edges of Bloch bands. But line 
solitons can also bifurcate from high-symmetry points inside Bloch bands, as has been reported 
numerically and experimentally before |13l 1141 . Here the high-symmetry points are points on 
the dispersion surface \x = fj,(k x ,k y ) where dfj,/dk x = d^/dky = 0, with k x ,k y being the Bloch 
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Figure 5: Eigenvalue comparison for offsite line solitons with a = 1 and lattice (15): (top) tan# = 
1; (bottom) tan# = 2. (Left) Numerically calculated A value in solid blue and the asymptotic 
approximation in dashed red. (right) The value for C left from taking the ratio of the numerical 
eigenvalues and the asymptotic expression as well as the predicted limiting value (red star) from 
equations ( 38 ) and ( 39 ) in section [5j 
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Potential 




Outer Equation 





Off Site 




Figure 6: (Color online) Line solitons in the asymmetric lattice (68) with a = 1 and tan# = — 1. 



(Top left) the potential (68). (Top right) the phase calculated from the "outer" integral equation 
( |30[ ); the asterisk marks the phase value 0o at the singularity kq = 2y/2. (Bottom) onsite and 
offsite line solitons for e = 0.35 as well as the predicted centers of the solitons (dashed). 
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wavenumbers in the x and y directions. For line solitons bifurcated from such interior points of 
Bloch bands, however, resonance with Bloch modes may occur. Such resonance excites Bloch-wave 
tails that are non-vanishing in the direction perpendicular to the line and hence makes the line 
"soliton" nonlocal |20j . To obtain true line solitons inside Bloch bands, this resonance must be 
absent, a requirement that poses strong restrictions on the angles of line solitons. Indeed, we will 
see below that inside Bloch bands, line solitons at only a couple of special angles are admissible. 
First we consider a concrete example with the lattice 

V(x, y) = 6 sin 2 x + 4 sin 2 y. (70) 

Inside the first Bloch band, fj, G [3.6080, 4.1565], of this lattice there is an X-symmetry point 

{k x0 , k y0 , Mo) = (0, 1, 3.9529), (71) 

whose Bloch mode is 7r-periodic in x and 27r-periodic in y. The dispersion surface near this X-point 
is saddle-shaped. For line solitons bifurcating from this X-point, /i = Mo + r/e 2 , 77 = ±1, and e <C 1. 
Taking e = 0.2, the level curves of the dispersion surface at fJ,(k x , k y ) = Mo + r]e 2 for rj = ±1 are 
displayed in Fig. [7] (left panels, solid lines). 

Suppose the angle of the line soliton with the x-axis is 9. Then at the m value of the line soliton, 
linear Bloch modes that are periodic along this line direction, with the period matching that of the 
X-point Bloch wave, are located in the wavenumber plane at the intersections of the parametrically 
defined line 

k x = kq sin 6 + k x o, (72a) 
k y = —ko cos 9 + kyQ, (72b) 

with the level curve 

fi(k x ,k y ) = mo + r]e 2 . (73) 

Existence of such intersections means that the line soliton is in resonance with Bloch modes. 

Inspection of the level curves in Fig. [7] shows that for almost all angles 9, the line (72), whose 



slope is — cot 6, always intersects those level curves for both r\ = ±1. The only exceptions are 
9 = 0, ±7r/4 for 77 = 1 and 9 = n/2 for rj = —1. In these cases, resonance is absent, thus true 



line solitons could exist. Recalling the conditions for envelope solitons below Eq. (12), we see that 



true line solitons with angles 9 = 0,±7r/4 may bifurcate from the X-point (71) under defocusing 
nonlinearity (<r = —1), and true line solitons with angle 9 = ir/2 may bifurcate from this X-point 
under focusing nonlinearity (a = 1). But where are the envelope positions of these bifurcating 
line solitons inside Bloch bands? To answer this question, it is again necessary to use exponential 
asymptotics. 

Strictly speaking, our exponential-asymptotics analysis in the previous sections was for line 
solitons bifurcating from edges of Bloch bands and residing outside of them (thus resonance with 
Bloch bands never occurs). But for line solitons with special angles inside Bloch bands, those 
which avoid resonance, our previous analysis can be applied without any additional work. The 
conclusion is that, for each of those special angles, two families of line solitons, one onsite and the 
other offsite, bifurcate out from this X-symmetry point; and their envelope locations are given by 



the formulae (58). Numerically we have confirmed these predictions. For instance, the onsite line 



solitons near the X-point (71), at angle 9 = 7r/4 for defocusing nonlinearity and at angle 9 = tt/2 
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for focusing nonlinearity, have been numerically obtained and displayed in Fig. [7] (right panels); 
and their envelope locations agree with those predicted by the formulae (58). Our analytical results 



are consistent with the numerical and experimental reports of line solitons inside Bloch bands in 
[HI Q31 QS] as well. 



H(k x , k y ) = no + e 2 
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(Color online) In-band line solitons in the lattice (70) near the X-symmetry point 
(0,1) at special angles for e = 0.2. Top: rj = 1 and a = — 1 (defocusing nonlinear- 
Mo + f?e 2 



ity). Bottom: r] = —1 and a = 1 (focusing nonlinearity). Left: level curves n(k x ,k y ) 



extended periodically (solid) and the line ( 72 ) (dashed) for values of angle 8 which admit true line 



solitons. Right: onsite line solitons at these special angles with \i = [iq + rje 



As explained above, inside Bloch bands, line solitons for most inclination angles are nonlocal 
(in the normal direction), in the sense that they comprise non- vanishing Bloch- wave tails far away 
from the central line. Then a natural question is, how can we determine the amplitudes of those 
non- vanishing Bloch- wave tails? It turns out that this question can also be treated in the framework 
of the exponential-asymptotics theory developed above. Specifically, we should first recognize that 



the kq value at the intersection of the line (72) and the level curve (73) also satisfies equation 



(22) for singularity locations (which can be easily checked). Thus this acq value of resonant Bloch 
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modes is also a real pole in the Fourier transform ip(x,y,K) of the nonlocal line-soliton solution 
i[)(x,y,W). Unlike the fourth-order poles (26) in our earlier analysis (see Eq. (35)), this resonant 
pole is simple (i.e., first-order). The residue of this resonant pole is also exponentially small in 
e, and its exact value can be computed from the same outer integral equation (30) by the same 
numerical algorithm (32) in the earlier text. From the residue of this resonant pole and the inverse 
Fourier transform ( 56 ), the amplitudes of non- vanishing Bloch-wave tails in the nonlocal line solitons 
will then be obtained. This problem of calculating amplitudes of non-vanishing Bloch-wave tails 
in nonlocal line solitons resembles that of calculating amplitudes of continuous-wave tails in the 
fifth-order Kortewegde Vries equation, the third-order NLS equation and other related equations 
[5j [211 E2J [231 Ull ESI [26] . But our treatment of directly solving the outer integral equation for 
the residues of resonant poles differs from those earlier one-dimensional works, which hinge on the 
recurrence relation for the coefficients of the Taylor expansion of the Fourier transform. 

The results from the above specific example hold qualitatively for all line solitons inside Bloch 
bands. Specifically, inside Bloch bands true line solitons exist only for very few (up to three) special 
angles due to the requirement of resonance suppression. For each of those special angles, there exist 
two line solitons, one onsite and the other offsite. 



12 Conclusion 

In this paper, we have presented what we believe is the first step toward a fully two-dimensional 
asymptotic theory for the bifurcation of solitons from infinitesimal continuous waves. For line soli- 
tons bifurcating from infinitesimal Bloch waves in a general two-dimensional periodic potential, an 
analytical theory utilizing exponential asymptotics is developed. For this two-dimensional problem, 
the previous approach of relying on a certain recurrence relation to solve the outer integral equation 
is no longer viable due to the presence of complex poles which are closer to the origin than the 
real poles. Instead, we solved this outer equation directly along the real line up to the real poles. 
This new approach not only overcomes the recurrence-relation limitation, but also simplifies the 
exponential-asymptotics process. 

Using this modified technique, we showed that from every edge of the Bloch bands, line solitons 
with any rational line slope bifurcate out; and for each rational slope, two line-soliton families 
exist. In addition, a countable set of multi-line-soliton bound states have been constructed analyti- 
cally. Furthermore, as a byproduct of this exponential-asymptotics theory, a certain linear-stability 
eigenvalue that bifurcates out of the origin at a band edge, is analytically obtained. Line solitons 
bifurcating from interior points of the Bloch bands were investigated as well, and it was shown that 
such solitons exist (inside Bloch bands) only for a couple of very special line angles due to reso- 
nance with the Bloch bands. These analytical predictions were compared with numerical results 
for symmetric as well as asymmetric potentials, and good qualitative and quantitative agreement 
was obtained. 

Throughout the analysis, the potential in Eq. ([I]) is taken to be in minimal-period orientation. 
If the potential is not aligned along that minimal-period orientation, all our results would still hold, 
except that the smallest of the real poles in Eq. ( [26] ) would not be those given by Eq. (28), but 
rather be a certain multiple of those numbers. 

In this article, we assumed that the potential in Eq. ([I]) has equal periods in x and y. If the 
x- and y-periods are not the same, say n in x and x n i n with x being the ratio between the two 
periods, then we can introduce the scaled y variable y = y/x- I n this scaled variable, the periods 
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of the potential are tt in both x and y, hence satisfying the assumption in this paper. Due to 
this y-scaling, the Laplacian V 2 in Eq. changes to d xx + X~ 2 dyy- For this scaled "Laplacian", 
all our analysis is still valid (except for very minor modifications), because our analysis does not 
rely on the equal coefficients in the Laplacian at all. Following this analysis, we shall find that 
from every edge of the Bloch bands, line solitons with slopes of any rational number divided by x 
bifurcate out; and for each of those slopes, two line-soliton families exist. We shall also find that 
from a high-symmetry point inside Bloch bands, line solitons at only a couple of special angles may 
bifurcate out. 

Finally, we point out that the analysis in this paper is developed for line solitons bifurcating 
from high-symmetry points of the Bloch bands where the Bloch mode is unique (see Assumption 
1). At certain points of the Bloch bands, however, the Bloch modes are not unique [27J. In such 
cases, nonlinear interactions between different Bloch modes would occur [27J. To treat line-soliton 
bifurcations from such Bloch-band points, the exponential asymptotics analysis of this article would 
need to be generalized. Such generalizations will be left for future studies. 
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Appendix: Line solitons in stripe lattices 

In this appendix, we consider line solitons in stripe lattices, i.e., lattices which vary only in 
one direction V(x,y) = V(x). We shall show that line solitons in such lattices only have complex 
poles and no real poles. In addition, the envelopes of these solitons can be arbitrarily located. Our 
conclusion will be that complex poles do not place restrictions on envelope locations. 

The bulk of the analysis remains the same as in the main text. We consider stationary solutions 
whose leading-order term is a Bloch- wave packet, 

ip(x, y, W) ~ eA(W - W )b(x). (74) 

The packet envelope A(W—Wq) has a sech- profile and varies only in the direction of W = e(x sin 8 — 
ycosO), and Wq = exo is the location of this envelope. 



The pole singularities of these line wavepackets are the values of kq where Eq. (22) holds, except 
that the operator Lq drops the y-dependence now. That is, 

L (f> + iK Li4> - nl(f) = 0, (75) 

and 

L = dl + fjL ~ V(x) , Li = 2V • [sin 6, - cos 9) . (76) 



Converting this equation into the eigenvalue problem (27), we find that all eigenvalues kq are now 
complex- valued and not real (excluding 0). This is illustrated in figure [8] (right panel). Here the 
lattice is chosen as V(x) = 6 sin 2 x and the line slope as tan(#) = 1. 

In this stripe lattice, we have numerically found line solitons. An example with a = 1 is shown 
in figure [| (left panel). Since the lattice is y-independent, so are the Bloch modes b{x). Then if 
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Figure 8: (Color online) A line soliton (left) and its singularity structure (right) for the stripe 
lattice V[x) = 6sin 2 x with tan# = 1 and a = 1. 



a line soliton exists with the envelope centered at a particular Wq value, from equation (74) it is 



clear that line solitons with the envelope centered at arbitrary Wq values would exist, because any 
variation in Wq may be compensated for by a shift in the y coordinate of W. 

From the above analysis, we conclude that in a stripe lattice, envelopes of line solitons can be 
arbitrarily positioned, and poles of those solitons are all complex- valued (off the real axis). 
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